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ABSTRACT 

Infrared surveys indicate that the dust content in debris disks gradually declines with stellar age. We sim- 
ulated the long-term collisional depletion of debris disks around solar-type (G2 V) stars with our collisional 
code. The numerical results were supplemented by, and interpreted through, a new analytic model. General 
scaling rules for the disk evolution are suggested. The timescale of the collisional evolution is inversely pro- 
portional to the initial disk mass and scales with radial distance as r 4 3 and with eccentricities of planetesimals 
as e~ 23 . Further, we show that at actual ages of debris disks between 10 Myr and 10 Gyr, the decay laws of 
the dust mass and the total disk mass are different. The reason is that the collisional lifetime of planetesimals is 
size-dependent. At any moment, there exists a transitional size, which separates larger objects that still retain 
the "primordial" size distribution set in the growth phase from smaller objects whose size distribution is already 
set by disruptive collisions. The dust mass and its decay rate evolve as that transition affects objects of ever- 
larger sizes. Under standard assumptions, the dust mass, fractional luminosity, and thermal fluxes all decrease 
as with £ = —0.3...— 0.4. Specific decay laws of the total disk mass and the dust mass, including the value of 
£, largely depend on a few model parameters, such as the critical fragmentation energy as a function of size, the 
primordial size distribution of largest planetesimals, as well as the characteristic eccentricity and inclination of 
their orbits. With standard material prescriptions and a distribution of disk masses and extents, a synthetic pop- 
ulation of disks generated with our analytic model agrees quite well with the observed Spitzer/MIPS statistics 
of 24 and 70 /im fluxes and colors versus age. 

Subject headings: circumstellar matter — planetary systems: formation. 



1. INTRODUCTION 

Since the IRAS discovery of the excess infrared emission 
around Vega by Aumann et al. (1984), subsequent infrared 
surveys with ISO, Spitzer and other instruments have shown 
the Vega phenomenon to be common for main-sequence stars. 
The observed excess is attributed to second-generation cir- 
cumstellar dust, produced in a collisional cascade from plan- 
etesimals and comets down to smallest grains that are blown 
away by the stellar radiation. While the bulk of such a de- 
bris disk's mass is hidden in invisible parent bodies, the ob- 
served luminosity is dominated by small particles at dust 
sizes. Hence the studies of dust emission offer a natural tool 
to gain insight into the properties of planetesimal populations 
as well as planets that may shape them and, ultimately, into 
the evolutionary history of circumstellar planetary systems. 

In recent years, various photometric surveys of hundreds 
of nearby stars have been conducted with the Spitzer Space 
Telescope. These are the GTO survey of FGK stars (Beich- 
man et al. 2005; Bryden et al. 2006; Beichman et al. 2006a), 
the FEPS Legacy project (Meyer et al. 2004; Kim et al. 2005), 
the A star GTO programs (Rieke et al. 2005; Su et al. 2006), 
the young cluster programs (Gorlova et al. 2006), and others. 
These observations were done mostly at 24 and 70 /im with 
the MIPS photometer, but also between 5 and 40 with the 
IRS spectrometer (Jura et al. 2004; Chen et al. 2006). Based 
on these studies, about 15% of mature solar-type (F0-K0) 
stars have been found to harbor cold debris disks at 70 /im. 
For cooler stars, the fraction drops to 0%-4% (Beichman et al. 
2006a). For earlier spectral types, the proportion increases to 



about 33% (Su et al. 2006). At 24 /im, the fraction of systems 
with detected excess stays similar for A stars, but appreciably 
decreases for FGK ones. Similar results in the sub-millimeter 
range are expected to become available soon from a survey 
with SCUBA and SCUBA2 on JCMT (Matthews et al. 2007). 
Preliminary SCUBA results for M dwarfs suggest, in particu- 
lar, that the proportion of debris disks might actually be higher 
than suggested by Spitzer (Lestrade et al. 2006). 

All authors point out a decay of the observed infrared ex- 
cesses with systems' age. However, the values reported for the 
slope of the decay, assuming a power-law dependence t~ a , 
span a wide range. Greaves & Wyatt (2003) suggest a < 0.5, 
Liu et al. (2004) give 0.5 < a < 1.0, Spangler et al. (2001) 
report a w 1.8, and Greaves (2005) and Moor et al. (2006) 
derive a « 1.0. Fits of the upper envelope of the distribution 
of luminosities over the age yield a « 1.0 as well (Rieke et al. 
2005). Besides, the dust fractional luminosity exhibits a large 
dispersion at any given age. 

In an attempt to gain theoretical understanding of the ob- 
served evolution, Dominik & Decin (2003) assumed that 
equally-sized "comets" produce dust through a cascade of 
subsequent collisions among ever-smaller objects. If this dust 
is removed by the same mechanism, the steady-state amount 
of dust in such a system is proportional to the number of 
comets. This results in an M/Mq m r/t dependence for the 
amount of dust and for the number of comets or the total mass 
of the disk. Under the assumption of a steady state, this re- 
sult is valid even for more complex systems with continuous 
size distributions from planetesimals to dust. Tenuous disks, 
where the lifetime of dust grains is not limited by collisions 
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but by transport processes like the Poynting-Robertson drag 
(Artymowicz 1997; Krivov et al. 2000; Wyatt 2005), follow 
M oc t~ 2 rather than MocT 1 . 

More recently, Wyatt et al. (2007a) lifted the most severe 
simplifying assumption of the Dominik-Decin model, that of 
equal-sized parent bodies, and included them into the colli- 
sional cascade. A debris disk they consider is no longer a 
two-component system "comets + dust". Instead, it is a popu- 
lation of solids with a continuous size distribution, from plan- 
etesimals down to dust. A key parameter of the description by 
Dominik & Decin (2003) is the collisional lifetime of comets, 
t. Wyatt et al. (2007a) replaced it with the lifetime of the 
largest planetesimals and worked out the dependencies on this 
parameter in great detail. Since the collisional timescale is in- 
versely proportional to the amount of material, r oc 1/M , 
the asymptotic disk mass becomes independent of its initial 
mass. Only dynamical quantities, i.e. the disk's radial po- 
sition and extent, the orbiting objects' eccentricities and in- 
clinations, and material properties, i.e. the critical specific 
energy and the disruption threshold, as well as the type of the 
central star determine the very-long-term evolution. 

Still, there are two important simplifications made in the 
model by Wyatt et al. (2007a): (i) the disk is assumed to be in 
collisional equilibrium at all sizes, from dust up to the largest 
planetesimals and (ii) the minimum specific energy needed to 
disrupt colliding objects is independent of their size. As a con- 
sequence of (i) and (ii), the size distribution of solids is a sin- 
gle power-law. To check how reasonable these assumptions 
are, realistic simulations of the disks with collisional codes 
are necessary (e.g., Thebault et al. 2003; Krivov et al. 2005, 
2006; Thebault & Augereau 2007). 

The aim of this paper is two-fold. First, we follow the 
evolution of debris disks with our elaborate numerical code 
(Krivov et al. 2005, 2006) to check the existing analytic mod- 
els and the assumptions (i) and (ii) they are based upon. Sec- 
ond, in order to make these numerical results easier to use, we 
develop a new analytic model for the evolution of disk mass 
and dust mass that relaxes both assumptions (i) and (ii) above. 

Section 2 summarizes the basic ideas and assumptions and 
describes our numerical model and the runs of the collisional 
code. In Section 3 the numerical results are presented and de- 
pendences of the collisional timescale on the disk mass, dis- 
tance to the star, and mean eccentricity of parent bodies are 
derived. In section 4, the analytic model for the evolution of 
disk mass and dust mass is developed. Section 5 analyzes the 
evolution of dust luminosities. In Section 6, we use the ana- 
lytic model to synthesize representative populations of debris 
disks and compare them with statistics of debris disks derived 
from the Spitzer surveys. A summary is given and conclu- 
sions are drawn in Section 7. 



2. NUMERICAL MODEL AND DESCRIPTION OF RUNS 

2.1. Basic Approach 

For all numerical runs in this paper, we use a C++-based 
collisional code (ACE, Analysis of Collisional Evolution). 
The code numerically solves the Boltzmann-Smoluchowski 
kinetic equation to evolve a disk of solids in a broad range 
of sizes (from sub-micrometers to about a hundred of kilo- 
meters), orbiting a primary in nearly-Keplerian orbits (gravity 
+ direct radiation pressure + drag forces) and experiencing 
disruptive collisions. Collisions are simulated with available 
material- and size-dependent scaling laws for fragmentation 
and dispersal in both strength and gravity regime. The cur- 
rent version implements a 3 -dimensional kinetic model, with 
masses, semi-major axes, and eccentricities as phase space 
variables. This approach automatically enables a study of the 
simultaneous evolution of mass, spatial, and velocity distribu- 
tion of particles. The code is fast enough to easily follow the 
evolution of a debris disk over Gyr timescales. A detailed de- 
scription of our approach, its numerical implementation, and 
astrophysical applications can be found in our previous papers 
(Krivov et al. 2000, 2005, 2006). 

2.2. Disruption Threshold and Critical Specific Energy 

An object is said to be disrupted in a collision, if the largest 
fragment is at most half as massive as the original object. If 
the impactor's relative velocity is so high that the ratio of im- 
pact energy and target mass exceeds the target's critical spe- 
cific energy, QJj,, the target (and the impactor) are disrupted. 
For small objects, this binding energy is dominated by ma- 
terial strength, and for larger objects, self-gravity takes over. 
Both regimes are usually described by a sum of two power 
laws (Krivov et al. 2005, Sect. 5.1, and references therein) 

where "s" and "g" stand for strength and gravity, respec- 
tively. The reported values of the coefficients A s and A g 
vary by more than one order of magnitude, and we took 
A s = A g = 5 x 10 6 erg/g, in agreement with the reference 
case for basalt given by Benz & Asphaug (1999). The expo- 
nents are 36 s = —0.3 and 36 g = 1.5 (corresponding to —0.1 
and 0.5 in the mass scaling). With these parameters, the two 
power-law components contribute equally at s m 316 m, and 
the lowest binding energy, the minimum is reached at 
s w 129 m. The influence of the choice of parameters on the 
resulting evolution is discussed in Sect. 4. 

For computational reasons, we refrained from including a 
treatment of cratering collisions in the runs. Note that these 
were not taken into account in previous studies of the long- 
term evolution of debris disks (e.g. Dominik & Decin 2003; 
Wyatt et al. 2007a) either. Thebault et al. (2003) and Thebault 
& Augereau (2007), who focused on shorter time spans, did 
include this non-disruptive type of collisions that lead to the 
continuous erosion of a target by small impacting projectiles. 
They found the effect to be dominant for particles in between 
100 fim and 1 cm for the case of the inner (3 Pictoris disk, 
while big, kilometer-sized objects in the gravity regime are 
mainly lost to disruptive collisions (see Table 4 in Thebault 
et al. 2003). However, including cratering can lower the life- 
time of large objects, especially when relative velocities are 
low and disruptive collisions are rare. Another caveat is that 
cratering collisions alter the shape of the wavy size distribu- 
tion at the lower end (Thebault & Augereau 2007), which af- 
fects the observable thermal fluxes. 
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2.3. Collisional Outcomes 

The distribution of sizes and the velocities of fragments in 
an individual (catastrophic) collision has been subject to stud- 
ies for decades. Laboratory work was done on high- velocity 
impacts on scales of millimeters and centimeters (e.g., Fuji- 
wara et al. 1977; Fujiwara 1986; Davis & Ryan 1990). Statis- 
tics on the mass distributions of observed asteroidal fami- 
lies and geometrical considerations (Paolicchi et al. 1996; 
Tanga et al. 1999; Tedesco et al. 2005) as well as gravito- 
hydrodynamic simulations of fragmentation and reaccumula- 
tion (Michel et al. 2002) cover the range of larger, kilometer- 
sized bodies. On small scales, the resulting size distributions 
show a strong dependence on impact velocity and seem to 
indicate a turn in the power law at fragment sizes around 
« 1 mm (or « 1% of the size of the used targets). The 
slope for objects above that size is steeper than the one for 
smaller objects (Davis & Ryan 1990). However, Thebault 
et al. (2003) found that the ratio of these two slopes and the 
size at which the slope changes influence simulation results 
only slightly. On kilometer and larger scales, the fragmenta- 
tion is influenced by gravitational reaccumulation of relatively 
small fragments onto bigger ones. Hence, bigger fragments 
(~ 100km) will be overabundant, and conversely, smaller 
fragments (<~ 1km) underabundant, compared to the under- 
lying distribution without gravity. The slopes of the size dis- 
tribution n(s) oc s~ p of kilometer-sized objects are poorly 
known. A wide range from p = 3.5 up to p = 9.0 has 
been reported. These deviations in the kilometer regime are 
most probably the severest caveat of the power-law approxi- 
mation, because they are independent of the actual material 
and caused only by gravity. Nevertheless, we assume that 
fragments follow a single power-law distribution nf rag (s) oc 
s~ 3 - 5 , expecting the influence on the final collisional steady 
state to be only moderate. 

2.4. Commons for All Runs 

All disk models presented here are set up around a star of 
solar mass and luminosity. Parameters of the central star af- 
fect the disk evolution in various ways. They determine the 
size limit for grain's blowout by radiation pressure and orbital 
velocities at a given distance, thereby altering impact veloc- 
ities and rates. For late-type stars, strong stellar winds may 
affect the dust dynamics (Augereau & Beust 2006; Strubbe 
& Chiang 2006). On the observational side, dust tempera- 
tures and brightnesses are influenced. Here, we focus on the 
scalings for a fixed spectral type (G2V), and not on scalings 
between different types. 

The disks themselves all share the same material proper- 
ties and shapes. We adopt the material, described by a bulk 
density p = 2.5 g/cm 3 , the radiation pressure efficiency of as- 
tronomical silicate (Laor & Draine 1993), and a critical frag- 
mentation energy as specified in Sect. 2.2. We switched off 
the Poynting-Robertson effect, which is unimportant for de- 
bris disks under study, as well as stellar wind drag, which 
plays only a minor role around G-type stars. The fragments 
produced in an individual collision are distributed according 
to a single power law, dA oc s~ 35 ds oc m _11 / 6 dm. A 
biggest fragment size is assumed to scale with specific im- 
pact energy to the power of 1.24 (for details, see Krivov et al. 
2006). The initial mass distribution is given by dA oc m~ g , 
with q = 1.87, a value that accounts for the modification of 
the classical Dohnanyi's (1969) q = 1.833 through the size 
dependence of material strength (see, e.g., Durda & Dermott 



TABLE 1 
Description of numerical 

RUNS. 



Run 


Distance [AU] 


Cmax 




Nominal runs 




ii-0.3 


7.5-15 


0.3 


i-0.3 


15-30 


0.3 


O-0.3 


30-60 


0.3 


oo-0.3 


60-120 


0.3 




Additional runs 




i-0.1 


15-30 


0.1 


i-0.2 


15-30 


0.2 


i-0.4 


15-30 


0.4 



1997). The particle masses range from 4.2 x 10~ 15 g, corre- 
sponding to a radius of 74 nm, to 4.2 x 10 21 g, corresponding 
to 74 km. The stepping between the 60 mass bins is logarith- 
mic with a factor of « 4 between neighboring bins. The initial 
radial profile of the particle density was given by a slope of the 
normal optical depth of —1.0. The initial total mass of each 
disk was set to 1 M e (earth mass). 

2.5. Specifics of Individual Runs 

We have made four "nominal" runs, each of which corre- 
sponds to a certain radial part of the disk between 7.5 and 120 
AU from the star (Table 1). In these runs we assumed ini- 
tial eccentricities of planetesimals to be uniformly distributed 
between e m i n = 0.0 and e max = 0.3, spanning three bins 
centered at 0.05, 0.15 and 0.25. In addition, three runs with 
altered maximum eccentricity of 0.1, 0.2, and 0.4 were made 
for the 15-30 AU ring. In all the runs, we assumed that or- 
bital inclinations are distributed between / m j n = e m ; n /2 and 
Imax = e max /2 in accordance with the energy equipartition 
relation I — e/2. 

3. NUMERICAL RESULTS AND SCALING LAWS 

3.1. Evolution of Disks of Different Masses 

A debris disk is said to be in a quasi-steady state or quasi- 
equilibrium if the amounts of particles with different sizes on 
different orbits, while changing with time (therefore "quasi"), 
stay constant relative to each other. For brevity, we will often 
omit "quasi" and use simply "steady state" or "equilibrium". 
To express the condition of a quasi-steady state formally, we 
can introduce a phase space, in which a dynamical state of 
each particle is characterized by a vector p. That vector may 
be composed, for instance, of coordinates and velocity com- 
ponents. Alternatively, p may represent the set of orbital el- 
ements of the object. Let n(p, s, f) be the number of objects 
with radii in [s, s + ds] at phase space "positions" [p, p + dp] 
that the disk contains at the time instant t. The assumption of 
a quasi-steady state can now be expressed as 

n(p,s,t) = h(p,s) f(t). (2) 

The total disk mass, 

M disk (t) = J J n(p,s,t)dpds, (3) 

can be rewritten as 

M diBk (t) = /(t) J Jh(p,s)dpds (4) 

or, setting /(0) = 1, 

M disk (t) - f(t)M , (5) 
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where Mq is the initial disk mass. As long as objects are both 
created and lost in two-particle collisions, their gain and loss 
rates are given by 



10° 



n(p,s,t) = J J J J [G(p,s,pi,si,p 2 ,s 2 



-i(pi, si, p 2 , s 2 )S(p - pi)S(s - si)] 

xn(pi,si) f(t) ra(p 2 ,s 2 ) f{t) 

x dpidsidp 2 ds 2 , (6) 

where the function G(p, s, pi, s\, p 2 , s 2 ) describes the gain 
in population p, s due to collisions between pi, si and p 2 , 
s 2 and the function L(pi, si, p 2 , s 2 ) accounts for the loss in 
population pi, si in collisions with p 2 , s 2 . The disk mass 
changes at a rate 



or 



M disk (t) = J f h(p,s,t)dpds 
Mdisk(t) = /(*) / / n(p,s)dpds. 



From Eqs. (6) and (7), we find that M disk (t) oc f 2 (t). 
Eq. (8) suggests M disk (t) oc f(t). Hence, /(i) oc 
Integration yields 

/- 1 



Using Eq. (5) we obtain 



M disk (t) 



t/T 



Mq 
1+t/T 



and 



(7) 

(8) 
while 

(9) 



(10) 



(11) 



where 1/C = Mq • r, i.e. the product of the initial mass 
and a characteristic time. This relation is invariant under the 
transformation (t, M d ; sk ) — * (t ■ x, M d i sk /a;), even if C is not 
constant. Therefore, the mass scale of a system under colli- 
sional evolution is inversely proportional to its timescale. For 
example, doubling the initial total mass halves the collisional 
lifetime of the system. All curves in the M disk (i) plots can be 
shifted along lines of equal t ■ M disk . 

Dominik & Decin (2003) used this approach and equated 
the characteristic time r with the collisional lifetime of their 
"comets". At the initial phase t <C t, Eq. (10) gives 

M disk (t) w M (1 - t/r) . (12) 

If the system is old enough so that t » r, the total mass 
will be just proportional to t _1 . Particles whose lifetimes are 
independent of the total mass are exempt from the asymptotic 
one-over-i behavior. Examples would be the /3-meteoroids 
that are blown out and small particles in disks tenuous enough 
for the Poynting-Robertson effect to be more efficient than 
collisions. The total mass of such particles is oc t~ 2 (Dominik 
& Decin 2003). 

As we have shown, for the systems that undergo a steady- 
state collisional evolution, the factor C in Eq. (11) (or r) 
should be constant. To check this, we evaluated C = 
— Mdisk/ M 2 isk for every two subsequent time steps of the nu- 
merical runs. The results are given in Fig. 1. 

Instead of being constant at later times, C decreases, 
roughly following a power law C oc ^-2/3... -4/5 j ne ex _ 
planation is simple: the systems did not reach an equilibrium 
where t » r or at least t « r during their lifetime. The 
evolution of the total mass in Fig. 2 demonstrates that as well. 
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FIG. 1. — The coefficient C from Eq. (11) as a function of time for four 
nominal runs. The total disk mass and time in the runs are scaled according to 
-^disk tx t _1 to compensate for the difference in dynamical timescale. Note 
that the near-constancy of C at the beginning of the evolution is an artefact 
of the double-logarithmic plotting. The double-linear inset shows that the 
decrease of C is fastest at earlier times. 




« 0.4 - 



10' 10° 
Time t [yr] 

FIG. 2. — The evolution of the total mass in the four nominal runs. Again, 
the plateau at the beginning of the evolution is an artefact of the logarithmic 
plotting of time. In fact, the mass decay is strongest at the very beginning 
(see inset and Eqs. (10), (12)). 

3.2. Dependence on Distance from the Star 

Rings of identical mass but at different distances have dif- 
ferent collisional timescales. The comparison in Fig. 1 shows 
that doubling the distance requires a 20-fold increase in disk 
mass to have the same timescale. This corresponds to a 
power-law dependence 



C oc r 



-4.3 



(13) 



In a thorough analytic approach based on a Dohnanyi-type 
collisional cascade, Wyatt et al. (2007a) came up with C oc 
r -i3/3 > wn ich is in good agreement with our numerical result. 
This index is made up of three contributions. First, the density 
in the rings drops with r~ 3 as their circumference, height, and 
width increase linearly. Second, the relative velocities have 
an r -1 / 2 dependence. Third, these impact velocities affect 
the minimum required mass for a projectile to be disruptive 
and thereby the total number of such projectiles. That gives 
another r 1_<? , where q is the slope in the appropriate mass 
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FIG. 3. — The influence of the average eccentricity of planetesimals on the 
timescale of disk's collisional evolution. Top: the evolution of the parameter 
C from Eq. (11) for four different runs (i-0.1, . . . , i-0.4). Bottom: four initial 
C values versus average eccentricity e = (e max + e m ; n )/2 (pluses) together 
with the C oc e 9 / 4 fit for those runs (line) and the same for runs with a 
narrower range of eccentricities, as described in Sect. 3.3 (crosses). 

distribution, e.g. q = 11/6 for the classical Dohnanyi case. 
See Sect. 4.3 for details. 

3.3. Dependence on Eccentricities of Parent Bodies 

The intrinsic collisional probability of planetesimals is 
nearly independent of their eccentricities, as long as they are 
not too high (see, e.g. Krivov et al. 2006). Nevertheless, ec- 
centricities determine impact velocities and, through that, the 
minimum size of a disruptive projectile. Therefore, higher 
planetesimal eccentricities imply a larger rate of catastrophic 
collisions and thus a faster collisional evolution. To quantify 
the dependence, we have made runs with maximum eccen- 
tricities of 0.1, 0.2, 0.3, and 0.4 (Table 1) and determined the 

values of C. The results suggest a power law C oc emax as 
shown in Fig. 3. 

This result comes as a surprise. Wyatt et al. (2007a) de- 
rive C oc e 5 / 3 . The same scaling is inherited by our analytic 
model, see Eq. (36) below. Since this discrepancy can be ei- 
ther due to an incompleteness of the analytic approach or to 
a non-linear relation between the maximum and the effective 
eccentricity, we tried to rule out the latter case by perform- 
ing additional runs with e confined to narrow bins of width 
0.1, centered at 0.05, 0.15, 0.25, and 0.35. These runs can 



be well described by the same power law, C oc e 9 / 4 (Fig. 3). 
Therefore, the analytic model fails to reproduce this particular 
dependence. Nevertheless, it correctly describes many others, 
as the next sections will show. 

4. ANALYTIC MODEL FOR EVOLUTION OF DISK MASS AND DUST 

MASS 

4.1. Size and Mass Distributions 

In what follows, we will analyze size or mass distributions 
of objects. Different authors use distributions of different 
physical quantities (number, cross section, mass) with differ- 
ent arguments (particles size or mass) and of different type 
(differential, cumulative, per size decade, etc.). A standard 
choice is to use a differential size distribution, n(s), that gives 
the number of particles per unit size interval: 



n(s) = J n(p,s)dp, 



(14) 



or a differential mass distribution, n(m), that gives the num- 
ber of particles per unit mass interval. Instead of n, it is often 
convenient to use the mass-per-size-decade distribution, 



dM, 



disk 



= ln(10) s m(s) n(s). 



(15) 



dlog 10 s 

In contrast to n(s), this quantity tells us directly, objects in 
which size range contribute the most to the mass of the sys- 
tem. Therefore, we will use it when plotting size or mass 
distributions. 

In the case of a power-law size distribution, n(s)ds oc 
s 2 ~ 3q ds is the number of objects with sizes [s, s + ds] and 
n(m)dm oc m~ q dm is the number of objects with masses 
[to, to + dm]. The mass per size decade is oc s 6 ~ 3q oc to 2-9 . 
When q < 2, the total mass is determined by large bodies, 
whereas the cross section is dominated by small particles as 
long as q > 5/3. 

4.2. Three-Slope Distribution 

The combination of material strength at smaller sizes and 
self-gravity at larger ones, with a turnover at around 100 m, 
causes the size distribution in a collisionally evolving system 
to strongly deviate from a single-slope power law, especially 
for object sizes of around 1 km. This is illustrated by Fig. 4 
that shows how a disk evolves from the first-guess power law 
to a more realistic size distribution. The speed of this evo- 
lution is determined by the collisional timescales of popu- 
lations of different-sized particles in the disk. Populations 
of smaller particles with sufficiently short lifetimes consist 
mostly of fragments of disruption of larger bodies. They will 
have reached collisional equilibrium with each other soon, ac- 
cording to their production rate by populations with longer 
lifetimes. Those latter populations of bigger particles will still 
be on their way to a steady state. As time goes by, more and 
more long-lived populations will undergo the transition from 
primordial to reprocessed material. 

As this transitional mass moves towards larger objects with 
time, the smaller particles follow to a new "intermediate 
steady state". The lower panel of Fig. 4 shows the develop- 
ment of the characteristic wavy shape in the size distribution 
(e.g., Campo Bagatin et al. 1994; Thebault et al. 2003; Krivov 
et al. 2006) at the small-size end near the blowout limit due to 
radiation pressure. Once established, this shape remains con- 
stant. Only the absolute level changes because this distribu- 
tion at smaller sizes acts as the trail of the distribution at larger 
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FIG. 4. — Results of the ii-0.3 ran. Top: Time evolution of mass in indi- 
vidual mass bins, from the largest bodies of 74 km in radius to the smallest, 
74 nm in radius. The mass ratio between adjacent bins is 4. Each solid line 
corresponds to one individual bin and gives the mass contained in that bin 
(see the right axis) as a function of time. The left axis can be used to find 
the line that corresponds to a given object size. The thick dashed curve cor- 
responds to fa 1 mm radius, i.e. to the largest solids still treated as dust. 
The thick dotted curve, which goes roughly through the turning points of the 
curves, is the transition size s t (t); see Eq. (37). Bottom: Size/mass distribu- 
tion at four specific instants of time shown in the top panel with vertical lines: 
initially, after 5 X fO 5 years when st has reached Sb, and after 5 X f0 7 and 
5 X If] 9 yr when significant dust depletion has already occurred. 

sizes. In the upper panel of Fig. 4, the number of smaller par- 
ticles is constant for some time and then goes down, as soon 
as the distribution in the gravity regime starts to deviate from 
its primordial one. 

These arguments suggest that an overall size distribution 
n(s) can be approximated by a combination of three power 
laws (Fig. 5). For particles large enough to be only barely af- 
fected by collisions at time t, we assume n to follow s 2 ~ 3,Jp . 
Here, q p is the "primordial" slope determined by the processes 
in which these planetesimals have formed. Small particles 
that are in quasi-steady state are separated from bigger pri- 
mordial objects by a transition zone which we characterize 
by a time-dependent size s t (t). To distinguish between the 
strength and gravity regimes, we introduce two more power 
laws and assume the mass distribution to follow n oc s 2 ~ 3,?B 
for gravity-dominated quasi-steady state and n oc s 2 ~ 3q " for 
strength-dominated quasi-steady state. The two regimes are 
separated by an object size Sb, which we will call breaking 
radius. Thus, the waviness is neglected, but the effect of a 
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FIG. 5. — Schematic plot of the three regimes in the mass distribution and its 
time evolution. The mass St divides second generation material in collisional 
equilibrium (s < St) from primordial material (s > St), while divides 
the material strength regime (s < s b ) from the gravity regime (s > s h ). 

size-dependent QJj is kept. 
The resulting size distribution is given by 



n(s) 

for s t < s < s max , 

n(s) = n mi 
for Sb < s < s t , and 

n(s) = n ma 



I Smax \ 
V S ) 



■St 



3« P -2 



St 



3(/ P 



3q g -2 



(16) 



(17) 



(18) 



for s m in < s < s b , where n max = n(s max ), with s max being 
the size of the largest planetesimals. From this distribution, 
two important quantities can be derived. One is the total disk 
mass, 



M disk 



nps ds, 



(19) 



and the other is dust mass (that determines the infrared lumi- 
nosity and therefore provides a link to observations), 



where s min < s d < s b . 



Sd 

/ n{s) t 



7rps 3 ds, 



(20) 



4.3. Collisional Lifetimes of Planetesimals 

As seen from Eqs. (16)-(20), the evolution of M d i s k and 
Mdust is controlled by n max (i) and s t (t). 

We start with n max and assume, according to Eqs. (2) 
and (9): 

^max(O) 



<(*) = 



1 + t/T„ 



(21) 



where T max is the collisional lifetime of these largest bodies. 
Equation (21) closely reproduces the disk evolution as soon 
as the whole system has reached the quasi-steady state at all 
sizes or, in other words, as soon as s t (t) has reached s max . 
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The second quantity that we need, s t (t), could easily be ob- 
tained by inverting the function r(s), the collisional lifetime 
of planetesimals of a given size s. To obtain r(s), we begin 
with the lifetime of the largest objects in a disk. Assuming 
that q > 5/3, Wyatt et al. (2007a, their Eq. 12) approximated 
it as 



7"max — 



4-7T 
CTtot 



•*i'/p 



r 5/2 dr 



(GM, 



,1/2 



X 77 V- (22) 

f(e,I)G(q,s) 

where e and / are the effective orbital eccentricities and in- 
clinations, CTtot is the initial cross sectional area of the disk 
material, G the gravitational constant, r the radial distance of 
the ring of parent bodies, and dr its width. The slope q in 
their single-power-law approach corresponds to the primor- 
dial slope q p in our nomenclature. The functions / and G are 
given by 



f(e,I) = 
G(q,s) = 



+ I 2 , 



(23) 



5-3q 



+ 2 



5/3 



+ 



9-4/3 
g-5/3 
q-l 



x c ( s f- 3c > - (^p) 



4-3 9 



3-3?' 



with 



X c (s) = 



GM* 



1/3 



(24) 



(25) 



While /(e, /) describes the dependence of the impact veloci- 
ties on eccentricities and inclinations, the functions G and X c 
characterize the disruption of planetesimals by smaller pro- 
jectiles. Namely, X c (s) is the minimum size ratio between 
the smallest disruptive projectile and the target, and G(q, s) is 
the number of disruptive projectiles. 

We need the lifetime of objects of an arbitrary size, t(s < 
s max ). To derive it, we can simply substitute s max by s in 
Eq. (22), obtaining 



4-7T 

r(s) = — 

CTtot 



3q P 



,5/2 dr 



(GM.) 



1/2 



X r. (26) 

fG{q p ,s) 

In order to replace the dependence on the initial cross sec- 
tional area of objects, otot, with their initial total mass, Mo, 
we need to derive both quantities from the initial size distri- 
bution in Eq. (16). The area is given by 

3q p -5 



CTtot 



c(0) 



7TS" 



3<?p - 

Since it is dominated by s„ 



CTtot 



c(0) 



1 



for q p > 5/3, we obtain 

■ a 3 /o \ 3q P -5 



3q p — 5 \ s 



The initial total disk mass is 



M = n max 



(0) 



3(6 - 3q p ) 



6— 3q p 



(27) 



(28) 



(29) 



For q p < 2, it is dominated by s max . However, since a primor- 
dial slope q p > 2 is not unrealistic (see Sect. 4.8) we refrain 
from using a further approximation. Then, the area and the 
mass are related through 



<7tot = M ■ 



1 



3(2 - gp) 
4(q p - 5/3) ' ' 

6— 3q p 



P 

^min 

S max 



3q P 



(30) 



Inserting Eq. (30) into Eq. (26) results in 

3q p — 5 



T(S): 



16irp 
3M 

g P - 5/3 
2-<? P 



,5/2 dr 



(GM*) 1/2 

6-3(j p 



f(e,I)G(q p ,sY 



(31) 



which gives the collisional lifetime of an object with radius ,s. 
Note that 



1 



g P 



1 - 



6— 3<3 P 



3 In- 



(32) 



for q p — > 2. 

If the mean impact velocities in the system are high enough 
to allow planetesimals of radius s to get disrupted in a colli- 
sion, i.e. X c (s) <C s m ax/s, G(q p , s) reduces to 



G(q p ,s) 



Qp ~ 5/3 
g P - 1 



X c (s) 



3— 3<j p 



(33) 



and t(s) to 



T(8): 



16irp 
3M ' 

g P - 1 
2-gp 



3g p 



r 2 dr • 



9P-1/2 



6 — 3<j p 



GM* , 



(34) 



Now, we take into account the dependence of on the ob- 
ject size s, as was done by O'Brien & Greenberg (2003). If 
we are only interested in the gravity regime, s > s D , Eq. (1) 
is simplified to 



Sb 



:s/, u 



(35) 



where QJj b is the critical specific energy at the breaking ra- 
dius, i.e. around the minimum of Q^(s). Assuming, further, 
that / oc e, we can write down the dependencies of the colli- 
sional lifetime, 

T(s) OC a to J • s 3<fe-5+3(9p-l)& g _ r 3/2+q p . dr . g-5/3 _ (3g) 

O'Brien & Greenberg (2003) yield the same size dependence 
on s in their Eq. (11). 

To find s t (t), the object size below which a steady state is 
reached, we assume that the populations move from their pri- 
mordial state to the quasi-steady state instantaneously when 
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the system age reaches their initial mean collisional lifetime, 
r(st) = t. Inverting that, the resulting mass of objects in tran- 
sition can be retrieved as a function of system age. Keeping 
the assumption X c <C s max /s, the relation is 

at^cxtVCaap-B+aCte-i)^) (37) 

for t > r(sb) = Tb. This transitional size is also plotted in 
Fig. 4. 

Pan & Sari (2005) followed a similar approach in their 
study of the Kuiper-belt size distribution. Describing the 
propagation of the shock wave through the target, they in- 
troduce a parameter (3 that varies between 3/2 (if all energy 
of a projectile goes to the shock wave) and 3 (if all its mo- 
mentum does). Their (3 equals l/b g in our nomenclature, and 
b g = 0.5 leads to /? = 2. Additionally, we have to replace 
their slope q with our 3q p — 2. Then, given their Eqs. (6), 
(7), and 7V >S oc s 3 ~ 3qp , we yield the same exponent as in 
our Eq. (37). Note that what Pan & Sari (2005) call "break- 
ing radius" is our "transition radius" s t , and their "radius of 
equilibrium" is our "breaking radius" Sb- 

4.4. Evolution of Disk Mass 

Now, we derive the full expression for the time-dependent 
total disk mass. Using the size distribution given by Eq. (18), 
fimax from Eq. (21) and expressing n max (0) through M with 
the aid of Eq. (29), we can perform the integration in Eq. (19). 
Then, the resulting time-dependent disk mass is 

r . « o. -i -1 

Mo 



M disk (i) : 



1 



+ 



6 — 3q p 



6-3q p 



3g p 



3g p 



\Smax ) 
\ ''max / 



6-3<? p 



6— 3<jp 



3q s -3q p 



2-9 P _ 2-q p 
2-g 3 



( s «) 



6— 3(3 a 



(38) 



for Tb < t < T max . To make Eq. (38) valid for earlier phases, 
i.e. for t < Tb, Sb should be replaced by St(t). The sizes 
involved are the maximum object size s max , the transition 
size between the primordial and reprocessed material s t , the 
breaking radius between the gravity and strength regime Sb- 
The lower limit in the size distribution, s m j n , is crucial for the 
dust emission and it is usually taken to be the radiation pres- 
sure blowout limit. As long as q p < 2, it is fairly unimportant 
for the mass budget. However, we are interested in q p > 2 
as well. Therefore, we can safely set s m ; n = only in the 
last line of Eq. (38), where it enters through s m j n /s max to the 
power of 6 — 3q s , with q s « 11/6 < 2. 

The relative importance of the terms in Eq. (38) is illus- 
trated in Fig. 6. A combination of the classic Dominik-Decin 
behavior in the first line of Eq. (38) together with the second 
line is a reasonably accurate approximation to Mdi s k(i) for 
most of the time. With the aid of Eq. (37), Eq. (38) trans- 
forms to 

-l 



M disk (i) f 



Mr, 



1+t/l 



6-3(j p 



Sb 



6-3<2 p 



^-9p 

9p-5/3+(g p -l)b g 
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FIG. 6. — The contributions of different terms in Eq. (38) (dotted and 
dashed lines) and their total (solid line). 
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FIG. 7. — Evolution of total masses with (scaled) time, obtained in four 
numerical runs and with the analytic model. 



1 



2 - , 



(39) 



for Tb < t < T max . At t -C T max , and assuming q p = 1.87, a 
further approximation is 



M disk (t) w M (1 - const • i a2 ) 



(40) 



The evolution of the disk mass, both from the numerical 
runs and from the analytic solution (38), is plotted in Fig. 7, 
showing a good agreement between analytics and numerics. 
A deviation is only seen around t = Tb where the transi- 
tion from primordial to reprocessed state sets on for gravity- 
dominated objects. The reason is that, to ease the analytic 
treatment, we neglect the smooth natural transition from ma- 
terial strength to self-gravity given by Eq. (1) and assume a 
sharp break between the two power laws instead. 

4.5. Evolution of Disk Mass at Latest Stages 

As soon as the age of the system has reached the collisional 
lifetime of the largest bodies, i.e. at t > T max , the solids of all 
sizes in the disk reach quasi-steady state, and the change in 
total mass will be dominated by \jt. At this latest phase, the 
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projectiles that can destroy objects of size s max no longer fol- 
low a size distribution with the primordial slope, 2 — 3q p . In- 
stead, they have the slope of a collisional cascade under grav- 
ity regime, 2 — 3q g . The slightly longer collisional lifetime 
can neither be expressed through Eq. (22) that uses the ini- 
tial cross section <7 tot nor through Eq. (31) that contains the 
initial disk mass Mq and slope q p . The correct way to eval- 
uate r max is to use the initial number density of biggest ob- 
jects, n max (0), and the slope q g . Expressing tr tot in Eq. (22) 
through n max with the help of Eq. (28) and replacing then q p 
with q g , we obtain 

12<7„ - 20 r 5 / 2 dr 



c(0) 



(GM«) 



1/2 



I 



• (41) 

f(e,I)G(q e , 

max ) 

Expressing now n max (0) through M by virtue of Eq. (29) 
yields 

_ IQirp r 5 / 2 dr 



3M 



(GM.) 



1/2 



5/3 



2 - 



9p 



6— 3(j p 



(42) 



where both slopes, g p and g g , appear (cf. Eqs. 22 and 31). 

4.6. Evolution of Mass in Dynamically "Cold" Disks 

All the treatment above applies to planetesimal belts where 
relative velocities are high enough for the biggest objects to 
be destroyed by mutual collisions. This might not be the case 
in dynamically "cold" disks with low eccentricities and incli- 
nations and/or very far from the star. 

Consider again the lifetime of objects r(s). As s increases, 
X c (s) (Eq. 25) increases too and at a certain point reaches 
s max /s. At this point, G (Eq. 24) becomes zero and r(s) 
(Eq. 31) goes to infinity. This means that, for a given impact 
velocity, objects above a certain critical size cannot be dis- 
rupted anymore. In systems with low relative velocities, that 
critical size may happen to be smaller than s max . This will 
affect the mass evolution. Specifically, when s t reaches that 
critical size, the overall mass decay ceases. 

To illustrate such effects, Fig. 8 shows the influence of the 
effective e and I on the evolution of the total mass of a disk 
of initially 1 M ffi at an effective distance of 10 AU, calculated 
with our analytic model. For colder disks, the curves start to 
flatten. This happens because the largest planetesimals (that 
dominate the total mass) stay intact, which slows down the 
mass loss. 

4.7. Evolution of Dust Mass 

The dust mass can be evaluated in a similar way as the disk 
mass. We use now Eqs. (18), (20), (21), (29), and (37). Ne- 
glecting the minimum mass s m j n only when it enters the for- 
mula through s m ; n /s max , we obtain 



M dust (t) = 



1 + t/r, 



Sb 



max 
2-9 P 



t 

Tb 

£d 
Sb 



9p-5/3+(g p -l)b g 



2-q s 



Sb 



-1 
(43) 
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FIG. 8. — Influence of the effective eccentricity assumed in the analytic 
model for a disk of 1 at r = 10 AU with a radial extent dr = 7.5 AU. 
The I = e/2 relation between eccentricity and inclination is assumed. 
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FIG. 9. — Similar to Fig. 7 but for dust masses, i.e. masses in particles with 
radii below 1 mm. 

for Tb < t < r max . Before that, i.e. at t < Tb, we have q s 
and b s instead of q g and b g , respectively. If the assumed pri- 
mordial slope, <7 p , equals the steady-state slope in the strength 
regime, q s , the dust mass stays constant, which is the case for 
the first part of the numerical integration. However, as soon 
as the transitional zone reaches objects large enough to be in- 
fluenced by self-gravity, Eq. (43) starts to work. It shows that 
the evolution of dust mass depends most strongly on the dif- 
ference between q p and q g . The dust mass decay, obtained 
both from the numerical runs and analytic solution (43), are 
shown in Fig. 9. For t > Tb, we roughly have Md us t oc 
with £ w -0.3. 

We finally note that Eq. (43) is valid as long as the colli- 
sional lifetime of the largest planetesimals is longer than the 
age of the system. When t > T max , t/Vb in that equation must 
be replaced by T max / T b . 

4.8. The Model Parameters 

Our analytic model contains several parameters that either 
differ from similar parameters in the numerical model (such 
as e) or are absent there (such as q s and q g ). To use the ana- 
lytic model, we have to specify them. We now describe how 
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FIG. 10. — Index £ of the power-law evolution of the dust mass, M d <x t«. 
The horizontal axis gives the dependence on the slope of the primordial mass 
distribution, q p , for values from cj g = 1.57 (bottom) to q s = 1.77 (top) 
for the slope in the gravity regime. The bold line is for q g = 1.67 ~ 5/3. 
Vertical lines indicate the mean value and error estimates for q p from Trujillo 
et al. (2001). 

this can be done, explaining, in particular, the choice of pa- 
rameters used to plot analytic curves in Figs. 6-9. 

Two important free parameters of the analytic model are 
q s and q g . We use the work of O'Brien & Greenberg (2003) 
who found the slope of the size distribution in a system in a 
collisional steady state. With the dependence of the critical 
specific energy on the object size given in Eq. (1), they give a 
power-law index 



q = 



1 1/6 + 6 
1 + 6 



(44) 



in their Eq. (24). With b = b s = —0.1 for the strength regime 
we have q = q s = 1.877. Similarly, with b = b g = 0.5 for the 
gravity regime, Eq. (44) can be used to derive q g w 5/3. It is 
these values that we used in Eq. (38) to produce Figs. 6-8 and 
in Eq. (43) to plot Fig. 9. 

In contrast to q s and q g , the primordial slope, q p , is a free 
parameter not only in the analytic model, but also in the nu- 
merical one. As stated in Sect. 2.4, in all "nominal" runs we 
assumed q p = 1.87, which corresponds to p p = 3q p — 2 = 
3.61 in the size scaling. In principle, q p describes the mass 
distribution at the onset of the collisional grinding of the disk 
and, therefore, represents a link to the planetesimal formation 
process. The outcome of the agglomeration phase is the input 
to the phase of disruptive collisions. The Kuiper belt is the 
only source for observational constraints to this parameter so 
far, and recent surveys suggest a value of p p = 4.0 ± 0.5 (e.g., 
Trujillo et al. 2001; Bernstein et al. 2004) or q p = 2.00+0.17. 
Simulations by Kenyon & Bromley (2004) yield p p = 4.0- 
4.5 or q p — 2.00-2.17. According to Eq. (43), where we have 
Mdust cx t*, and together with q g « 1.67, this would change 
the dust mass evolution from Md us t oc t~ - 32 for q p = 1.87 to 
M dust oc t^ 0A0 for q p = 2.00. Fig. 10 shows the rather mod- 
erate dependence of the index £ on the two mass distribution 
slopes, <7 g and q p . 

While the dust size limit, s d > has little influence on the mass 
budget, the breaking size, Sb, the maximum size, s max , and 
the ratio of the two are relevant to the evolution as they define 
the lifetime of the largest bodies r max relative to t\>. What 
is more, the ratio Sb/s max determines the rate of the mass 



decay in Eq. (39). From Sect. 2.2 we know the location of 
the breaking radius to be 316 m for the material properties 
assumed, and the upper size limit of all the runs was set to 
Smax = 74 km. 

Another parameter in the analytic model is the collisional 
lifetime of objects of breaking radius, Tb = r(mb). Eq. (31) 
expresses it through other parameters critical for the efficiency 
of collisions: the radial distance to the star r, the disk ra- 
dial extension dr, and the effective eccentricity e and inclina- 
tion /. We choose to fix both the effective distance and the 
disk extension to be r = 4/3dr = 10 AU when reproduc- 
ing analytically the results of the ii-0.3 run, 20 AU for i-0.3, 
40 AU for O-0.3, and 80 AU for oo-0.3. Further, the incli- 
nation can be coupled to eccentricity by assuming the equi- 
librium condition I = e/2. Thus, only e remains as a free 
parameter. The best fit to, e.g., the ii-0.3 run is achieved if 
we assume e w 0.075 in the analytic model, which is approx- 
imately one quarter of e max = 0.3. With these choices, we 
find r(sb) ~ 4 x 10 5 years. 

Alternatively, Tb can be directly retrieved from the break 
in the evolution of the dust mass (see Fig. 9). This method 
gives r(sb) « 5 x 10 5 years, which is approximately 4/3 
times the value calculated with Eq. (31). This discrepancy is 
probably a result of the particle-in-a-box assumptions made 
by Wyatt et al. (2007a) in derivation of Eq. (22). We prefer 
this empirical scaling and thus applied the factor of 4/3 to all 
analytically estimated timescales in this paper. 

5. EVOLUTION OF DISK LUMINOSITY 

5.1. Fractional Luminosity for a Given Age 

Following Wyatt et al. (2007a), we define the fractional lu- 
minosity of dust as 

/ d = a tot /(4^r 2 ), (45) 

which assumes that dust grains are black bodies, absorbing 
and re-emitting all the radiation they intercept. Wyatt et al. 
(2007a, their Eq. 20) found that there is a maximum possi- 
ble fractional luminosity / max for a given age, whose value 
is independent of the initial disk mass, but depends on other 
model parameters such as the distance r of the disk center 
from the star, its width dr, size of the largest planetesimals 
D c , critical fragmentation energy QJ,, orbital eccentricity of 
planetesimals e (with their inclination being J = e/2), as well 
as the stellar mass M* and luminosity L. t . 

We now wish to explore /d(t) and check whether it has an 
upper limit in the framework of our analytic model. To this 
end, we used Eq. (45) and calculated <7 to t with the aid of our 
Eq. (43) for the dust mass. We assumed a solar-type star with 
M* = L* = 1 and probed disks with Mdisk = 1, 3, 10, and 
30M e ; r = 3,10,30, and 100 AU; dr/r = 1/8,1/4,1/2, 
and 1; e = 0.05,0.10,0.15, and 0.20. The results are pre- 
sented in Fig. 1 1 (thick lines). As a standard case, we adopted 
M disk = 10M®, r = 30 AU, dr/r = 1/2, and e = 0.10. It is 
shown with a thick solid curve in each of the panels. 

In the same Fig. 11, we have overplotted with thin lines 
the dust luminosity / d computed with Eqs. (14), (19), and 
(20) of Wyatt et al. (2007a), for comparison. In that calcula- 
tion, we assumed Q{j = 300 J/kg (constant in their model), 
D c = 60 km, and the same values of those parameters that are 
common in their and our model (M», L*, r, dr/r, and e). 

Analysis of Fig. 1 1 allows us to make a number of con- 
clusions. First, as expected, our model yields more gently 
sloping curves than that by Wyatt et al. As discussed above, 
the 1 jt law will be asymptotically reached in our model, too, 




FIG. 11. — Fractional luminosity of dust around a solar-like star as a function of age. Thick lines: our analytic model; thin lines: / d of Wyatt et al. (2007a). 
Different panels demonstrate dependence on different parameters: M disk (top left), r (top right), dr/r (bottom left), and e (bottom right). A standard case with 
M, = L* = 1, M disk = IOMq, r = 30 AU, dr/r = 1/2, and e = 0.10 is shown with solid lines (common in all panels). 



but this does rarely happen at ages t < 10 Gyr. Only the 
first signs of the curves' steepening appear at Gyr ages, and 
that only for the cases when the collisional evolution is faster 
(higher masses, closer-in or more confined dust rings, higher 
eccentricities). As a consequence of the slope difference be- 
tween the two models, our model places more stringent upper 
limits of /d at earlier ages, and conversely, it allows the Gyr- 
old systems to have a somewhat higher /<j than the model by 
Wyatt et al. does. 

Next, the dependence of / max on the initial disk mass, 
which cancels out in their model, is retained in our nominal 
runs (the top left panel). In fact, the maximum possible /<j is 
then determined by the maximum initial disk mass that still 
appears physically plausible in the framework of theories of 
planetesimal accretion and planet formation. 

Another point to mention is that, whereas the dependence 
on the disk width (bottom right) and planetesimal eccentrici- 
ties is relatively weak and monotonic, the dependence on the 
disk location (top right) is rather strong and more intricate. 
That the dependence is strong is the consequence of Eq. (13) 
that predicts the timescales to very sensitively depend on the 
distance from the star, and of Eq. (45) that contains a "dilution 
factor" r 2 . At the beginning of the evolution, the innermost 
ring is always the brightest because the dilution factor r 2 in 
Eq. (45) is the smallest. At the end of the evolution, the op- 



posite is true: the outermost ring will become the brightest, 
because its collisional evolution is the slowest and it retains 
more mass than inner disks. Therefore, all four curves inter- 
sect each other at a certain point; the 30 and 100 AU curve do 
that after 10 Gyr, i.e. outside that right edge of the plot. After 
that, all the curves go parallel to each other in the "Dominik- 
Decin regime", following a 1/t law. Note that inner rings 
reach the 1/t regime more quickly: already at 10 AU it is 
established in around 100 Myr for an initial mass of 10 Mq. 

Although the existence of a "maximum fractional luminos- 
ity for a given age", as suggested by Wyatt et al. (2007a), no 
longer holds in our model as a robust mathematical statement, 
in practice our model still suggests that fd(t) cannot exceed 
a certain limit, unless the model parameters take extreme val- 
ues, incompatible with our understanding of the planetesimal 
disks. For instance, we do expect /d < 10~ 4 at t — 10 Gyr, 
provided that the initial disk did not contain more than 30 
earth masses of solids and that the mean orbital eccentric- 
ity of planetesimals is not lower than 0. 1 (corresponding to 
the mean inclination larger than 3°). Therefore, plots such 
as Fig. 1 1 can be used to check whether or not /<j observed 
for a certain system with a known age is compatible with a 
"smooth", unperturbed collisional evolutionary scenario. In 
case it is not, it will be an indication that other mechanisms 
(delayed stirring, recent giant break-ups, non-collisional dust 
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FIG. 12. — Flux ratio versus time for (top) 24 and (bottom) 70 fim. 

production etc.) should be thought of to explain the observa- 
tions. 

5.2. 24 and 70 Micron Fluxes from Partial Rings 

In order to produce directly observable quantities from the 
derived dust masses, we now concentrate on dust luminosi- 
ties at particular infrared wavelengths. We calculated the 
dust temperature and the thermal emission integrated over 
the whole disk with a more accurate, yet sufficiently sim- 
ple, model, assuming that the absorption/emission efficiency 
is constant up to wavelengths of 2ir times the size of the par- 
ticles, s, and proportional to s _1 beyond that (Backman & 
Paresce 1993). Then we computed the spectral flux densities 
of dust emission and of the stellar radiation F* at a certain 
wavelength, as well as their ratio Fd/F*. As the size distri- 
bution in the dust regime quickly reaches its steady state, the 
luminosity F^ is directly proportional to the dust mass. There- 
fore, the same initial constancy and subsequent decay with 
£ = -0.3 ... - 0.4 apply. 

Fig. 12 shows the evolution of the excess emission at the 
Spitzer/MIPS wavelengths 24 and 70 /im, obtained from the 
four nominal runs. Since all disks have the same initial total 
mass (1 M®), the disks closer to the star are brighter and start 
to decay earlier. The difference between the excesses at 24 
and 70 [im, a measure of the disks' effective temperature, is 
varying with radial distance as well. Thus, the convergence 
of just the 70 fim fluxes at later times is only coincidental. 




Time [yr] 



FIG. 13. — Time evolution of the infrared excess of extended disks with 
different initial radial distributions (labels indicate the radial slope of the sur- 
face mass density; the thicker lines, the flatter the profiles) at 24 fim (dashed 
lines) and 70 (im (solid lines). The total mass is lMj in each case. 

It is a result of the radial dependence of temperature and the 
collisional timescale. 

5.3. Fluxes from Extended Disks 

Since resolved debris disks suggest that the parent body 
reservoir in the disks is usually confined to a toroidal region (a 
planetesimal belt), or is made up of several such tori, it seems 
appropriate to simply combine individual rings without taking 
into account possible interactions between particles that be- 
long to different rings. Thus, we summed up the fluxes from 
the four main runs. Different radial distributions in the whole 
disk were simulated by "weighting" the individual rings: 

4 

F d = ^F d;j (r,/r r, (46) 
i=i 

where rj are the central distances of the rings and values of 0, 
1, 2, and 3 were used for the slope 7. As the reference runs 
were made for rings of one earth mass each with volumes pro- 
portional to r 3 , the corresponding volume density in the ex- 
tended disk is proportional to r 7 " 3 , while the pole-on surface 
density and normal geometrical optical depth follow oc r 7 ~ 2 . 
The distance ro normalizes the total mass to 1M®. Therefore, 
by changing the slope, the mass is only shifted between inner 
and outer regions. 

In Fig. 13 the effect on the 24 and 70 fim fluxes is shown. 
If the weights are assigned in favor of more distant debris 
rings, the resulting fluxes are naturally reduced. The same 
is true for the speed of the decay because the timescales get 
longer. The evolution of the fluxes at the two Spitzer/MIPS 
wavelengths 24 and 70 yum differs significantly. At 24 fim the 
decay starts earlier and reaches its maximum speed earlier be- 
cause shorter-lived inner regions make the main contribution. 

The models contain a sufficient number of parameters, vari- 
ation of which would affect the curves in Fig. 13 in different 
ways. As stated earlier, varying the total mass changes the 
timescale according to r oc M^ isk . Hence, the curves can be 
shifted along the lines of equal t ■ Md; s k, i.e. along the top 
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left - bottom right diagonal. As seen from Fig. 13, variation 
of the radial distribution changes both the absolute level and 
the tilt of the curves. Besides, it affects the disk colors, i.e. 
the separation of the 24 and 70 /im curves in Fig. 13. In ad- 
dition, the dynamical timescales, and therefore the tilt of the 
curves, are affected by eccentricities and inclinations of the 
parent bodies that may reflect the presence of planetary per- 
turbers in the disk (see Sect. 3.3). Altogether, these degrees 
of freedom would allow one to reproduce a broad set of ob- 
servational data. 

6. COMPARISON WITH OBSERVATIONAL DATA 
6.1. Spitzer Data 

The advent of the Spitzer Space Observatory has brought 
a tremendous increase in the number of main-sequence stars 
surveyed for the existence of cold dust emission (see Werner 
et al. (2006) for a recent compilation). 

The wealth of data from these debris disk surveys allows us 
to confront our models with actual observations. To this end, 
we searched the literature for published flux ratios at 24 and/or 
70 /im (two of the three MIPS bands) around G-type main- 
sequence stars. To qualify as a main-sequence star we applied 
a lower limit to the stellar age of 10 Myr. Sources with stel- 
lar age estimates younger than this are likely stars with gas- 
dominated, protoplanetary disks; these were not taken into 
account. 

The bulk of the data taken in the framework of the 
Legacy program "Formation and Evolution of Planetary Sys- 
tems" (FEPS) (Meyer et al. 2004, 2006) is public since De- 
cember 2006. The FEPS archive contains images, spec- 
tra, photometry tables and Kurucz photosphere models 
and is available at http://data.spitzer.caltech.edu/popular/feps/ 
20061223_enhanced_vl/. Age estimates have been published 
for 46 FEPS G stars (Kim et al. 2005; Stauffer et al. 2005; 
Silverstone et al. 2006). 

The large Guaranteed Time Observer (GTO) survey of FGK 
stars contains another 64 stars, where ages are available (Be- 
ichman et al. 2005, 2006a; Bryden et al. 2006). Data for ten 
more G stars are listed in Chen et al. (2005a,b). In total, 120 
G-type main-sequence stars with flux ratios at 24 and/or 70 
fxm have been compiled from the literature for comparison 
with model flux ratios. 

6.2. Population Synthesis 

Based on the analytic prescription presented in Sect. 4 and 
motivated by the Wyatt et al. (2007b) work, we now build a 
synthetic set of debris disks around G2 stars. We generate 
a set of ring-like disks of width dr located at distances r G 
[r m i„,r max ], with masses M d i s k G [M min ,M max ], and ages 
between 10 Myr and 10 Gyr. The probability to have a disk of 
initial mass M at radius r was assumed to follow M " 1 r -0 - 8 , 
where M^ 1 corresponds to a log-normal distribution of initial 
disk masses and the r~ - 8 dependence was proposed by Wyatt 
et al. (2007b). As described in Sect. 5.2, the temperatures and 
the resulting thermal fluxes are calculated using the modified 
black-body formulas by Backman & Paresce (1993), assum- 
ing the emitting grains to have s = 1 /im, in agreement with 
the size distribution shown in Fig. 4. The other parameters are 
taken to be: q p = 2.00, q g = 1.67, q s = 1.877, dr/r = 0.5, 
21 = e = 0.15, Qfj(l m) = QJj(l km) = 5 x 10 6 erg/g, 
&d = —0.12, b g = 0.47, roughly corresponding to basalt in 
Benz & Asphaug (1999). Due to the small observational sam- 
ple, our aim was not to perform a multi-parameter fit to the 
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FIG. 14. — Flux ratios versus time for 24 [im (top) and 70 [im (bottom). 
The synthesized population (small dots) is compared to the observed one 
(big dots). Individually labeled is the possibly transient system HD 72905, 
see text. 



observations, but rather to cover the range of observed flux 
densities, which is defined by the limits of the distributions, 
not by their slopes. 

Varying disk locations and masses easily reproduces the ob- 
served distribution of fluxes at 24 /xm and 70 /im (Fig. 14). 
The synthetic population shown corresponds to r m ; n w 
20 AU, r max w 120 AU and M min < 0.01 M e , M max w 
30 Mq. Here, the radial range is needed to cover the range of 
colors, i.e. the ratios between the excess emissions at the two 
wavelengths. The mass range is needed to cover the observed 
range of excess, especially for younger disks at 70 /zm. 

Analyses of Spitzer detections might indicate a statistically 
significant increase of both 24 and 70 pm fluxes at ages be- 
tween a few tens of Myr to a few hundreds of Myr. (e.g., 
J. M. Carpenter et al., in prep.), which can only be marginally 
seen in our sample (Fig. 14). It is hypothesized that this fea- 
ture is caused either by an increased dust production due to 
delayed stirring by growing planets or by events similar to 
the late heavy bombardment in the solar system. Such effect 
could only be studied with an improved version of our ana- 
lytic model or with the numerical one. 

The distribution of disk colors is more difficult to repro- 
duce. Fig. 15 shows a significant abundance of fainter but 
warmer disks in an area that is not covered by the synthetic 
population. One explanation would be that the upper mass 
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FIG. 15. — Relation between fluxes at 24 and 70 fim versus time. The 
synthesized population (small dots) is compared to the observed one (big 
dots and triangles). The average photospheric uncertainty for both filters is 
marked by dashed lines in the upper panel. Excesses below those limits in 
either of the two filters are marked by triangles in the lower panel. In addition, 
the upper panel shows lines of equal dust mass and the lower panel gives the 
ring radii corresponding to the colors. 

limit is a function of radial distance, and that the innermost 
disks tend to be less massive and less luminous, from the very 
beginning. In addition, the lower panel of that Fig. 15 shows 
a trend towards higher effective temperatures for higher ages, 
which is difficult to understand. Indeed, as long as faint close- 
in disks are observed around older stars, one would expect 
ever brighter disks, and therefore more numerous detections 
of disks at the same distances around younger stars. Further- 
more, the trend in question contradicts to the results by Najita 
& Williams (2005), who found no significant correlation be- 
tween the disk radii and ages. Most likely, the discrepancy 
is only caused by uncertainties of the measured excesses at 
24 /im. Bryden et al. (2006) report that the average photomet- 
ric accuracy in that filter band is only as good as 1(724 = 6% 
due to stellar photosphere fitting errors and flat-field uncer- 
tainties. Therefore, excesses below those 6% of the pho- 
tospheric emission cannot be considered as significant. For 
70 /im, Bryden et al. (2006) state la 7a w 15%. Both limits 



are shown in the upper panel of Fig. 15. 

In Figs. 14 and 15, there is one particular system directly 
labeled. That system, HD 72905, was observed to show sig- 
nificant excess emission not only at 24 and 70 /urn, but also in 
the spectral ranges 8-13 and 30-34 iim of the Spitzer/IRS in- 
strument (Beichman et al. 2006b). The presence of two dusty 
regions was suggested: one exozodiacal at 0.03-0.43 AU and 
one around 14 AU. From the excess at 8-13 /im, Wyatt et al. 
(2007a) inferred the dust population in HD 72905 to be tran- 
sient because the observed fractional luminosity is above the 
maximum expected for a system of 300-400 Myr. As long as 
only 24 and 70 /im are considered, the HD 72905 dust does 
not seem particularly hot or bright, although it is among the 
hotter disks. 

At this point, it is interesting to compare our results to those 
of Wyatt et al. (2007b). Both analytic approaches aim at ex- 
plaining and reproducing the observations. Our model is dif- 
ferent from theirs in that we take into account the size de- 
pendence of the critical specific energy as well as the tran- 
sition from a "primordial" size distribution of planetesimals 
to the one set up by a collisional cascade. The amount of 
dust in their model is determined, from the very beginning, 
by the rather long collisional timescales of objects of tens of 
kilometers, so that the collisional evolution is much slower. 
This can be seen from the equations: in the denomina- 
tor of Eq. (10) causes the mass to stay almost at the initial 
level for a long time, before the system reaches the i -1 de- 
cay. In our model, although the mass decay is asymptotically 
slower (t> with £ w —0.3 ... — 0.4), it sets up very quickly, 
namely on collisional timescales of objects with minimum 
binding energy (st, ~ 100 meters). Therefore, we would 
expect the model by Wyatt et al. (2007b) to show signifi- 
cantly larger excesses at ages considered, if all other parame- 
ters were comparable. This, however, is not the case. Wyatt 
et al. (2007b) assumed a much weaker material in their colli- 
sional prescription. Their QJj = 300 J/kg at an object radius 
of 30 km (D c ~ 60 km) is by more than two orders of mag- 
nitude below the values we use in Eq. (1). As r a 1 
in Eq. (34), their collisional timescales are shorter and their 
evolution faster, too. Besides the material strength, the dif- 
ference in the assumed effective eccentricities — e = 0.05 in 
their model against e max /2 = 0.15 in ours — causes another 
factor of roughly 10 in the collisional timescales, according 
to Sect. 3.3. All the differences listed happen to nearly com- 
pensate each other. As a net result, the excesses predicted by 
Wyatt's et al. and our models are comparable with each other 
(see also Fig. 11), being in reasonable agreement with the ob- 
served ones. 

7. SUMMARY AND CONCLUSIONS 

We investigated the long-term evolution of debris disks 
around solar- type (G2V) stars. Firstly, we performed nu- 
merical simulations with our collisional code. Secondly, 
the numerical results were supplemented by, and interpreted 
through, a new analytic model. The latter is similar to, and 
builds up on, the model developed earlier by Wyatt et al. 
(2007a), but extends it in several important directions. It natu- 
rally includes the transition from the "primordial" size distri- 
bution of left-over planetesimals, set up at their agglomeration 
phase, to the size distribution established by the collisional 
cascade. Further, it lifts the assumption that the critical spe- 
cific energy needed for disruption is constant across the full 
range of sizes, from dust to the largest planetesimals. With 
these improvements, a good agreement between the numerics 
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and analytics is achieved. 
We draw the following conclusions: 

1. The timescale of the collisional evolution is inversely 
proportional to the initial disk mass. For example, halv- 
ing the total mass doubles all collisional timescales. 
This rule is valid for systems where collisions are the 
only loss mechanism of particles and only as long as /3- 
meteoroids are unimportant for the collisional budget. 

2. Numerics and analytics consistently yield a r oc r 4 3 
dependence of the timescale of the collisional evolution 
on the radial distance. 

3. Numerical simulations show that the collisional 
timescale varies with the average eccentricity of dust 
parent bodies as r oc e -23 . The analytic approach sug- 
gests a somewhat weaker dependence, r oc e -5 / 3 . 

4. An evolving three-slope size distribution is proposed to 
approximate the numerical results. The biggest objects 
are still distributed primordially, with a slope q p . The 
objects below a certain transitional size are already re- 
processed by collisions and thus have a quasi-steady- 
state size distribution, determined by their self-gravity 
(for intermediate-sized objects, slope q g ) or by ma- 
terial strength (for smallest objects, slope q s ). That 
transitional size corresponds to the largest objects for 
which the collisional lifetime is still shorter than the 
age of the system. The transitional size increases with 
time, meaning that ever-larger planetesimals get in- 
volved into the collisional cascade. 

5. At actual ages of debris disks, ~10 Myr to ~10 Gyr, 
the decay of the dust mass and the total disk mass fol- 
low different laws. The reason is that, in all conceivable 
debris disks, the largest planetesimals have longer col- 
lisional lifetimes than the system's age, and therefore 
did not have enough time to reach collisional equilib- 
rium. If the system were let to evolve for sufficiently 
long time, both dust mass and disk mass would start to 
follow i -1 . However, this requires time spans of much 
longer than 10 Gyr. 

6. The loss rate of the dust mass, and the decay rate of 
fractional luminosity, primarily depend on the differ- 
ence between the slope q p of the "primordial" size dis- 
tribution of largest planetesimals and the slope q g of 
the size distribution of somewhat smaller, yet gravity- 
dominated, planetesimals that already underwent suffi- 
cient collisional evolution. With "standard" values of 
q p and q g , the dust mass and the thermal fluxes follow 
approximately Ifi with £ = —0.3 ... — 0.4. 

7. Specific decay laws of the total disk mass and the dust 
mass largely depend on a few model parameters. Most 
important are: the critical fragmentation energy as 
a function of size, the slope of the "primordial" size dis- 
tribution of planetesimals q p and their maximum size 
s max , and the characteristic eccentricity e and inclina- 
tion I of planetesimals. 

8. The property that the maximum possible dust lumi- 
nosity for a given age does not depend on the initial 
disk mass, established by Wyatt et al. (2007a), is only 



valid in cases of very rapid collisional evolution, i.e. in 
closer-in or dynamically very hot disks. For most of the 
systems at ages < 10 Gyr, an increase of the initial disk 
mass leads to an increase of the dust luminosity, unless 
that initial mass is assigned extreme values, incompati- 
ble with our understanding of planetesimal disks. 

9. Assuming standard material prescriptions and disk 
masses and extents, a synthetic population of disks gen- 
erated with our analytic model generally agrees with the 
observed statistics of 24 and 70 fluxes versus age. 
Similarly, the synthetic [24] -[70] colors are consistent 
with the observed disk colors. 

As every model, our numerical model makes a number of 
general simplifying assumptions; the analytic one imposes 
further simplifications: 

• The collisional evolution is assumed to be smooth and 
unperturbed. Singular episodes like the aftermath of gi- 
ant break-ups or special periods of the dynamical evo- 
lution such as the late heavy bombardment are not in- 
cluded. 

• Effects of possible perturbing planets are taken into 
account only indirectly: through the eccentricities of 
planetesimals (dynamical excitation) and confinement 
of planetesimal belts (truncation of disks). Further ef- 
fects such as resonant trapping or ejection of material 
by planets are neglected. 

• We only consider disruptive collisions. This is a reason- 
able approximation for disks that are sufficiently "hot" 
dynamically. However, cratering collisions become im- 
portant when the relative velocities are insufficient for 
disruption to occur. 

• Neither dilute disks under the regime of Poynting- 
Robertson drag nor very dense disks with collisional 
timescales shorter than orbital timescales and with 
avalanches (Grigorieva et al. 2007) are covered by the 
present work. 

• Explaining the initial conditions or deriving them from 
the dynamical history of the systems at early stages 
of planetesimal and planetary accretion was out of the 
scope of this paper. Correlations between disk masses, 
disk radii, and the presence of planets, for example, 
were not considered, although they might alter the scal- 
ings we found here. 

Despite these limitations, our models reproduce, in essen- 
tial part, the observed evolution of dust in debris disks. We 
hope that they may serve as a starting point for in-depth stud- 
ies that will certainly be undertaken in the future, motivated 
by questions that remain unanswered, as well as by new data 
expected from ongoing and planned observational programs. 

We wish to thank Jean-Francois Lestrade, Philippe 
Thebault, and Mark Wyatt for fruitful discussions and Amaya 
Moro-Martm for useful review comments. This research has 
been partly funded by the Deutsche Forschungsgemeinschaft 
(DFG), project Kr 2164/5-1. 
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